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ABSTRACT 

Large Scale Shocks are responsible for the heating of the ICM and can be important 
sources of Cosmic Rays (CR) in the Universe. However the occurrence and properties of 
these shocks are still poorly constrained from both a theoretical and an observational side. 

In this work we analyze the properties of Large Scale Shocks in a (103Mpc/h) 3 cosmo- 
logical volume simulated with the public 1.0.1 release of the ENZO code. Different methods 
to identify and characterize shocks in post processing are discussed together with their uncer- 
tainties. Re-ionization affects the properties of shocks in simulations, and we propose a fitting 
procedure to model accurately the effect of re-ionization in non-radiative simulations, with 
a post-processing procedure. We investigate the properties of shocks in our simulations by 
means of a procedure which uses jumps in the velocity variables across the cells in the simu- 
lations and this allows us to have a viable description of shocks also in relatively under-dense 
cosmic regions. In particular we derive the distributions of the number of shocks and of the 
energy dissipated at these shocks as a function of their Mach number, and discuss the evolu- 
tion of these distributions with cosmological time and across different cosmic environments 
(clusters, outskirts, filaments, voids). 

In line with previous numerical studies (e.g. Ryu et al.2003, Pfrommer et al.2006) rela- 
tively weak shocks are found to dominate the process of energy dissipation in the simulated 
cosmic volume, although we find a larger ratio between weak and strong shocks with respect 
to previous studies. The bulk of energy is dissipated at shocks with Mach number M w 2 
and the fraction of strong shocks decreases with increasing the density of the cosmic environ- 
ments, in agreement with semi-analytical studies in the case of galaxy clusters. 

We estimate the rate of injection of CR at Large Scale Shocks by adopting injection 
efficiencies taken from previous numerical calculations. The bulk of the energy is dissipated 
in galaxy clusters and in filaments and the flux dissipated in the form of CR within the whole 
simulated volume at the present epoch is « 0.2 of the thermal energy dissipated at shocks; 
this fraction is smaller within galaxy clusters. 

Finally we discuss the properties of shocks in galaxy clusters in relation with their dy- 
namical state. In these regions the bulk of the energy is dissipated at weak shocks, with Mach 
number w 1.5, although slightly stronger shocks are found in the external regions of merging 
clusters. 

Key words: galaxy: clusters, general - methods: numerical - intergalactic medium - large- 
scale structure of Universe 



1 INTRODUCTION 

Galaxy clusters store up to several 10 63 ergs in the 
form of hot baryonic matter, with typical temperatures 
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of several keV. This thermal energy is mostly the by- 
product of shock-heating processes occurred during 
the formation of cosmic structures. However detect- 
ing shocks in Large Scale Structures (LSS) is still ob- 
servationally challenging since they usually develop in 
the external regions of galaxy clusters, where the X- 
ray emission is faint. In a few cases, however, internal 
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shocks driven by the merging events have been dis- 
covered with typical Mach numbers « 1.5 — 3 (e.g. 
Markevitch et al. 1999; Markevitch et al. 2002; Bel- 
sole et al.2004; Markevitch & Vikhlinin 2007). Stud- 
ies concerning the convolution of simulated cosmolog- 
ical shocks into simulated X-ray images (e.g. Gardini 
et al. 2004; Rasia et al. 2006) have clearly shown that 
the apparent lack of shocks in the ICM may be due to 
the concourse of projections and mass-weighting ef- 
fects along the line of sight, which tend to bias towards 
lower values of the overall temperature of the ICM and 
to smooth temperature gradients. 

Shocks are important not only to understand the 
heating of the ICM but also because they may be 
efficient accelerators of supra-thermal particles (e.g. 
Sarazin 1999; Takizawa & Naito 2000; Blasi 2001). 
Non thermal activity in galaxy clusters is proved by ra- 
dio observations which show synchrotron emission in a 
fraction of merging clusters: Radio Halos , at the clus- 
ter center, and Radio Relics, elongated sources at the 
cluster periphery (e.g., Feretti 2005). Several mecha- 
nisms related to cluster mergers and to the accretion of 
matter can act as sources of non thermal components 
in galaxy clusters. Large scale radio emission in the 
form of giant Radio Halos may be powered by parti- 
cle re-acceleration by MHD turbulence injected in the 
ICM during energetic merging events (Brunetti et al. 
2001; Petrosian 2001; Brunetti et al. 2004,08). Strong 
shocks may accelerate supra- thermal electrons from 
the thermal pool and explain the origin of Radio Relics 
(Ensslin et al.1998), while high energy electrons accel- 
erated at these shocks can produce X-rays and gamma- 
rays via Inverse Compton scattering off CMB pho- 
tons (e.g. Sarazin 1999; Loeb & Waxmann 2000; Blasi 
2001; Miniati 2003). Relativistic hadrons accelerated 
at shocks can be advected in galaxy clusters and effi- 
ciently accumulated (Volk, Aharonian, & Breitschw- 
erdt 1996; Berezinsky, Blasi, & Ptuskin 1997) produc- 
ing an important non-thermal component which could 
be directly sampled by future gamma ray observa- 
tions (e.g., Blasi, Gabici & Brunetti 2007). Secondary 
particles injected in the ICM via proton-proton colli- 
sions may also produce detectable synchrotron radia- 
tion (e.g. Blasi & Colafrancesco 1999; Dolag & Enssil 
2000) and may be eventually re-accelerated by MHD 
turbulence yielding an efficient picture to explain Ra- 
dio Haloes (Brunetti & Blasi 2005). 

The energetics associated with the population of 
cosmic ray particles (CR) accelerated at shocks de- 
pends on the Mach number of these shocks (e.g. Kang 
& Jones 2002). The Mach number distribution of cos- 
mological shocks is thus important to understand CR 
in galaxy clusters. Semi-analytical studies pointed out 
that shocks that form during cluster mergers are weak, 
M ~ 1.5, being driven by sub-clumps crossing the 
main clusters at the free-fall velocity (Gabici & Blasi 
2003, Berrington & Dermer 2003). These approaches 
however are limited as they treat cluster mergers as 
binary encounters between ideally virialised spheri- 
cal systems. Therefore cosmological numerical simu- 
lations represent a necessary avenue to address this is- 
sue in more detail. First attempts to characterize shock 



waves in cosmological numerical simulations were 
produced by Miniati et al.(2000), by employing a set 
of eulerian simulations and a shock detecting scheme 
based on jumps in the temperature variable. Later 
works adopted more refined shock-detecting schemes 
and were more focused onto the distribution of en- 
ergy dissipated at shocks (Keshet et al.2003; Ryu et 
al.2003, Hallman et al.2003, Pfrommer et al.20063). 
Ryu et al. and Pfrommer et al. basically found that the 
bulk of shocks in the universe is made of relatively 
weak shocks, but also allow to constrain the popula- 
tion of stronger shocks forming in the external regions 
of galaxy clusters, were structures are not completely 
virialised. In these environments, strong shocks are fre- 
quent and may provide the bulk of the acceleration of 
CR in large scale structures (Ryu et al. 2003; Pfrom- 
mer et al. 2006). For this reason numerical simulations 
are important to study the non-thermal components in 
galaxy clusters (e.g. Miniati et al.2001, Miniati 2003; 
Keshet et al.2003; Pfrommer 2008). However, the iden- 
tification and characterization of shocks, as well as the 
calculation of the energy injected in the form of CR at 
these shocks, is difficult because of the complex dy- 
namics of large scale structures and because of the se- 
vere limitations in terms of physics and numerical res- 
olution that affect present cosmological simulations. 

In this paper we follow the approach of the sem- 
inal paper by Ryu et al. (2003), studying the shock 
wave patterns in LSS and characterizing them in a 
post-processing phase. This allows us to briefly dis- 
cuss the effect of different shock detection techniques 
and their dependence on the variation of underlying 
physical processes (e.g. re-ionization) in the simula- 
tion. We use the cosmological code ENZO (Bryan & 
Norman 1997), which treats gas dynamics with an Eu- 
lerian scheme and allows us to follow the assembly of 
LSS with sufficiently good spatial resolution. 

The outline of the paper is the following. In Sectfl] 
we provide a brief introduction to the ENZO Code, 
in Sect[3] we present our cluster sample and the main 
properties of cosmological structures in our simu- 
lations, and in Sect[4] we discuss the effect of re- 
ionization on the thermal properties of simulated cos- 
mic structures. In Sectj5]we provide the different meth- 
ods to characterize shocks in post processing and in 
Sectj6] we discuss their main source of uncertainties in 
the cosmological framework. In SectJ7]we present our 
results about the main shocks properties and about the 
injection of CR. The main conclusions of this work are 
given in SectJU In the Appendix we the effect of the 
resolution and of as on our results. 



2 NUMERICAL CODE - ENZO. 

A precise description of the behavior of the bayonic 
gas is crucial for the goals of the present work. In par- 
ticular, the numerical code adopted for our simulations 

1 After this paper was submitted, also Skillman et al.2008 provided results 
on shock waves in ENZO AMR simulations, which are broadly consistent 
with the above works in literature. 
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Table 1. Main characteristics of the simulations. 



Volume 


Resolution 


physics 


ID 


(145Mpc) 3 


125fcpc 


adiab. 


AD125 


(145Mpc) 3 


250fcpc 


adiab. 


AD250 


(145Mpc) 3 


500fcpc 


adiab. 


AD500 


(145Mpc) 3 


SOOkpc 


adiab. 


AD800 


(80Mpe) 3 


125kpc 


cool. + reion. 


C0125 


(80Mpe) 3 


250fcpc 


cool., reion. and crS = 0.74 


S8250 


(80Mpe) 3 


250fcpc 


cool.+reion. 


CO250 



must support an accurate treatment of the dynamics of 
high supersonic flows and the formation and propaga- 
tion of strong shock waves during the process of cos- 
mological structures formation. The ENZO code sup- 
ports such description. ENZO is an adaptive mesh re- 
finement (AMR) cosmological hybrid originally writ- 
ten by Greg Bryan and Michael Norman (Bryan & 
Norman 1997, 1998; Norman & Bryan 1999; Bryan, 
Abel, & Norman 2001, Norman et al.2007). It cou- 
ples an N-body particle-mesh solver with an adap- 
tive mesh method for ideal fluidynamics (Berger & 
Colella, 1989). ENZO adopts an Eulerian hydrody- 
namical solver based on the the Piecewise Parabolic 
Method (PPM, Woodward & Colella, 1984). The PPM 
algorithm belongs to a class of schemes in which an 
accurate representation of flow discontinuities is made 
possible by building into the numerical method the cal- 
culation of the propagation and interaction of non- 
linear waves. It is a higher order extension of Go- 
dunov's shock capturing method (Godunov 1959). It 
is at least second-order accurate in space (up to the 
fourth-order, in the case of smooth flows and small 
time-steps) and second-order accurate in time. The 
PPM method describes shocks with high accuracy and 
has no need of artificial viscosity, leading to an opti- 
mal treatment of energy conversion processes, to the 
minimization of errors due to the finite size of the cells 
of the grid and to a spatial resolution close to the nomi- 
nal one. In the cosmological framework, the basic PPM 
technique has been modified to include the gravita- 
tional interaction and the expansion of the universe. For 
a detailed review of the ENZO code and a comparison 
with other numerical approaches, we refer to O'Shea 
et al.(2004) and O'Shea et al.(2005). 

In this work, in order to keep our study of LSS 
shocks as simple as possible, we use ENZO with a 
fixed spatial resolution without the application of AMR 
techniques. 



3 COSMOLOGICAL SIMULATIONS AND TESTS. 
3.1 General Properties 

In our simulations we have assumed a "Concordance" 
model, with density parameters Qq = 1.0, Qbm 
0.01 1, Q DM = 0.226, n A = 0.73, Hubble parame- 
ter h = 0.71, a power spectrum produced according 
to the Eisenstein & Hu (1999) fitting formulas with 
a primordial spectrum normalization erg = 0.94, and 
an initial redshift of z = 50. In order to have a large 
cluster statistics we simulated a total volume equiva- 
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Figure 1. Cumulative Mass Function for the total matter of all halos in the 
simulations, with poissonian errors. Press & Schechter (dashed) and Sheth 
& Tormen (dotted) mass functions are reported for comparison. 
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Figure 2. Total Mass versus Emission-Weighted Temperature for the sim- 
ulated galaxy clusters of our AD125 run (diamonds). Points for 15 clusters 
of the C0125 run are also plotted (asterisks). Best fit relations for these 
samples are drawn (solid line = AD125, dashed = C0125), together with a 
comparison the results of Borgani et al.2004 (yellow line). 

lent to (U5Mpc) 3 w (103Mpc/hf at the fixed nu- 
merical resolution of 125 kpc. This total volume was 
obtained by combining together six (independent) sim- 
ulated boxes of 80 Mpc per side. 

A list of all simulations used in our study with their 
main properties is listed in Tab. 1 . The goal of this study 
is to investigate cosmological shocks with the most 
simple physical and numerical treatments. Cosmolog- 
ical shock waves are supposed to be mainly driven 
by the assembly of cosmic structures, and therefore 
gravity should be the principal ingredient to model. 
Therefore we made massive use of "adiabatic" simu- 
lations at various resolutions (AD 125, AD250, AD250, 
AD800): these simulations contain only "adiabatic" 
physics, i.e. they do not have radiative cooling, UV 
photo-ionization at early epochs, thermal conduction 
and magnetic fields. These simulations are the start- 
ing point to investigate the effects on the properties 
of shocks driven by the adoption of a more complex 
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physical modeling. In particular, the re-ionization pro- 
cess has the important effect of increasing the tem- 
perature (and the sound speed) of cosmic baryons in 
the low temperature regions, and thus this is the first 
additional ingredient to take into account. Therefore 
we re-simulated at two spatial resolutions one of our 
six 80 Mpc boxes with the Haardt & Madau (1996) 
re-ionization model plus radiative cooling {CO 125, 
CO250) and used the outputs of these simulations to 
extract a recipe to mimic the effect of re-ionization in 
post-processing in adiabatic simulations. Finally, we 
perform simulations with a different <t 8 parameter, in 
order to study how this parameter can affect our results 
(S8250 simulation). 

We discus s the effect of re-ionization and cooling 
in Sectiond4l& l6.3l and the effect of the numerical res- 
olution and of erg in Appendix (Sects.A and B). 



3.2 Properties of the Simulated Galaxy Clusters 

The aim of this Section is to present the sample of 
galaxy clusters obtained from our simulations and to 
briefly discuss their main properties. 

A cluster reconstruction procedure, based on total 
over-density criteria (e.g. Gheller, Moscardini & Pan- 
tano 1998), has been applied to the outputs of the var- 
ious simulations at different cosmological times, pro- 
viding a population of synthetic galaxy clusters which 
can be followed during time. The overall AD 125 sim- 
ulation at z = consists of 113 galaxy clusters with 
total virial masses ^ 3 • 1O 13 M //i. 

The cumulative mass function of all the halos in 
our sample is reported in Fig fTJ and shows an overall 
good agreement with the Sheth & Tormen (1999) mass 
function for M < 3- 10 14 M o //i. The relevant deficit of 
halos with M ^ 3- 10 14 M //i is likely due to the artifi- 



cial cut-off in the over-density power spectrum at long 
wavelengths in our 80 Mpc simulations (Bagla & Ray 
2005; Bagla & Prasad 2006). Massive galaxy clusters 
are expected to be the most important regions in which 
kinetic energy is dissipated by shocks (in thermalisa- 
tion and CR acceleration). Because a mass deficit in 
cluster halos may introduce a deficit in the energy pro- 
cessed by shocks in the total simulated volume, this 
should be taken into account when we compare our re- 
sults with previous studies (e.g. Ryu et al. 2003; Pfrom- 
mer et al.2006; Kang et al.2007). 

In Figf2] we report the scaling between the total 
mass and the emission-weighted temperature within 
rsoo derived for our population of clusters Q Points are 
slightly above the self similar scaling found by Borgani 
et al.(2004) although they are consistent within a ~ 15 
per cent scatter; this is true also for the additional 15 
halos in the C0125 data set (with radiative cooling and 
re-ionization). 

Overall FigsfT] and [2] show that the basic statisti- 
cal and physical properties of our galaxy clusters are 
in line with those from other cosmological numerical 
simulations. 



4 PHASE DIAGRAMS AND RE-IONIZING 
BACKGROUND. 

Shocks and compressions driven by the gravitational 
force are the only sources that increase the thermal en- 
ergy of cosmic baryons in our adiabatic simulations. 
Baryons far away from collapsing regions have the 

2 7"500 is the radius encompassing a mean total over-density of 500 times 
the cosmic mean density, and roughly corresponds to 0.5 the virial radius of 
galaxy clusters in a ACDM cosmology. 
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Figure 5. Time evolution of the temperature distribution from z = 2 to z = (see panels for color-coding), for the C0125 simulation (Left) and for the 
corresponding adiabatic simulation AD125 with the post-processing treatment of re-ionization (Right). 
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Figure 4. Time evolution of the best-fitting relation for the minimum tem- 
perature of the COI25 run, from 2 = 2 down to z = 0. Best fit parameters 
for these curves are reported in Tab.2. 



Table 2. Fit parameters for the minimum temperature for the outputs of the 
C0125 ran. 

cl c2 



Redshift 



LogT 



: 2.0 
1.0 

: 0.8 

: 0.6 

: 0.4 

: 0.2 

: 0.0 



3.2383 ± 0.0032 
2.9388 ±0.0010 
2.8846 ± 0.0009 
2.8288 ±0.0011 
2.7628 ±0.0019 
2.6889 ± 0.0016 
2.5904 ±0.0016 



0.3198 ±0.0061 
0.5056 ± 0.0020 
0.5358 ± 0.0016 
0.5773 ± 0.0022 
0.5437 ± 0.0038 
0.5517 ±0.0033 
0.5711 ±0.0032 



0.0749 
0.0335 
0.0361 
0.0267 
0.0330 
0.0409 
0.0469 



± 0.0025 
± 0.0008 
± 0.0007 
± 0.0008 
± 0.0015 
± 0.0013 
± 0.0001 



lowest temperature that can be potentially affected by 
the process of re-ionization which occurred between 
z ~ 30 and z ~ 6 (e.g. Fan, Carilli & Keating 
2006). This process heats up the medium, increas- 
ing the speed of sound in the lowest temperature re- 
gions and this affects our estimate of the Mach num- 
ber of shocks. Therefore a study of cosmological shock 
waves must deal with the influence of a re-ionizing 



background (Haiman & Holder 2003; Loeb & Barkana 
2005; Mellema et al.2006). 



The re-ionization scheme available in ENZO is 
linked with a treatment of radiative cooling, which is 
computed by assuming an optically thin gas of pri- 
mordial composition, in collisional ionization equilib- 
rium, following Katz, Weinberg & Hernquist (1996). 
The time-dependent UV background is introduced ac- 
cording to Haardt & Madau (1996) and it models the 
effect of a population of quasars that re-ionizes the uni- 
verse at z « 6. The implementation of run-time re- 
ionization is more expensive in terms of memory usage 
compared to non-radiative simulations, and we thus 
applied it only in two re-simulated data-sets {CO 12 5, 
CO250). The effect of a re-ionizing background can 
also be modeled in the post-processing phase by in- 
creasing the temperature of each cell in the simulation. 
This has been done in Ryu et al.(2003), where a con- 
stant value of T = W 4 K was imposed to the simu- 
lated volume at z = 0. This may correctly reflect the 
complete re-ionization inside halos at present epoch 
(Haiman & Holder 2003), however it may overestimate 
the mean temperature of baryons far away of the most 
massive cosmic structures. Figure [3] shows the phase 
diagram in one (80Mpc) 3 simulation from the AD125 
data set, and in its re-simulation with cooling and re- 
ionization, CO 125. Re-ionization efficiently removes 
the coldest phase of the baryons and a forbidden region 
in the log T — log p space forms (where T and p are gas 
temperature and density, respectively), which actually 
traces the lower bounds of the Warm Hot Intergalac- 
tic Medium (Katz, Weinberg & Hernquist 1996; Cen 
& Ostriker 1999; Valageas, Schaeffer & Silk 2002; Re- 
gan, Haehnelt & Viel 2007). This lower bound is evolv- 
ing with time, as shown if Fig]4]where we report the fit 
to the value of the 15 per cent percentile in the distri- 
bution of temperature in the cells for different density 
bins. We also checked that variations in the percentile 
(up to ~ 50 per cent) does not significantly affect the 
results. By restricting to baryon densities in the range 
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1(T 32 < p < 3 • W~ 30 gr/cm 3 we obtain best fits with 
a second order polynomial : 



log(- 



Cl log(^)+c 2 (log(^)) 2 , 
Po Po 

-32, 



(1) 



where po = 10~ gr cm . The best fit parameters 
for each redshift are reported in Table 2. At moderate 
redshift (z ^ 1) all curves can also be approximated 
with a simple power law, T m i n oc p°' e (consistent with 
Valageas et al.2002), and with a normalization decreas- 
ing with time. In particular, at z = the minimum tem- 
perature of baryons is given by: 



T mm (K) = 450 (^) a60 , 
Po 



(2) 



which is indeed consistent with Eq.2 1 in Valageas 
et al. (2002). By using the CO250 data set we also veri- 
fied that these fits do not change with spatial resolution. 

We use the fitting formulas in Tab. 2 to increase the 
temperature of baryons in our adiabatic simulations in 
the post-processing analysis. In Figf5]we show the evo- 
lution with time of the temperature distribution func- 
tion, in case of adiabatic sub-sample of the AD 125 
simulation with post-processing treatment and in the 
case of C0125 simulation with run-time re-ionization. 
The agreement between the two distributions demon- 
strates the validity of our post-processing approach 
that will be applied in the following to the full set of 
our adiabatic simulations (AD 125, AD250, AD500 and 
AD800). 
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Figure 6. Maps of reconstructed Mach numbers using the VJ method based 
on jumps of 1 cells (top left), 2 cells (top right), 4 cells (bottom left) and 8 
cells (bottom right). The image is 20Mpc per side, the width along the line 
of sight is 125fepc. 
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5 SHOCK-DETECTING METHODS 
5.1 Basic Relations 

The passage of a shock in a simulated volume leaves its 
imprint as a jump in all the thermodynamical variables. 
Under the simple assumption that the pre-shocked 
medium is at rest and in thermal and pressure equi- 
librium, the pre-shock and post-shock values for any 
of the hydrodynamical variables (density, temperature 
and entropy) is uniquely related to the Mach number, 
M = v s /c s , v s being the shock speed in the region 
and c s the sound speed ahead of the shock itself. The 
Rankine-Hugoniot jump conditions contain all the in- 
formation needed to evaluate M; if the adiabatic index 
is set to 7 = 5/3 one has the well known relations (e.g. 
Landau & Lifshitz 1966): 



P2 

Pi 

T2 

Ti 

S2 
Si 



AM 2 
M 2 + 3' 

(5M 2 - 1)(M 2 +3) 



and 

(5M 2 



16M 2 



1)(M 2 + 3) M 2 + 3 



16M 2 



AM 2 



2/3 



(3) 
(4) 

(5) 
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Figure 7. Distribution of temperature fluctuations (dot-dash lines) and ve- 
locity fluctuations (solid lines) for non-shocked cells in the simulation, at 
four different over-density regimes (see labels within the panel). Data are 
taken from one box of side SOMpc of our AD125 run. 
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with indices 1,2 referring to pre and post-shock 
quantities, respectively, and where the entropy S is 

S = T/p 2 / 3 . 

The Mach number can be obtained from thejumps 
in one of the hydro dynamical variables (Eqs[3f{5]) or 
from a combination of them. It is well known that the 
value of the density jump saturates at pil P\ = 4 in 
the case of relatively large Mach numbers (Eqj3]) and 
thus strong shocks cannot be constrained from density 
jumps. For this reason temperature and entropy jumps 
provide the most effective tools to measure M. 

Eqsf3H5]describe shock discontinuity in the case of 
an ideal shock. In practice measuring the Mach number 
of the shocks in simulations is more problematic. Mat- 
ter falling in the potential wells may have chaotic mo- 
tions and the temperature distribution is usually patchy 
due to the continuous accretion and mixing of cold 
clumps and filaments into hot halos. All these complex 
behaviors establish velocity, temperature and density 
discontinuities across the cells in the simulated box. In 
a post-processing analysis this is expected to modify 
irreparably the strength of the jumps in the thermody- 
namical variables in the shocked cells with respect to 
that expected in the ideal case (Eqsf3]43]). Consequently 
the estimate of the Mach number from these equations 
is subject to unavoidable uncertainties (see Sec. 6. 2). 



5.2 The Temperature Jumps Method 

The analysis of jumps in temperature is commonly 
adopted to measure the strength of shocks in Eule- 
rian cosmological simulations (e.g., Miniati et al. 2001 ; 
Ryu et al. 2003). 

Cells hosting a possible shock pattern are prelimi- 
nary tagged by means of two conditions: 

• VT • VS > 0; 

• V • v < 0; 

An additional condition on the strength of the tem- 
perature gradient across cells is also customary re- 
quested, e.g. 

• | AlogT |> 0.11; 

(specifically | AlogT |^ 0.11 filters out shocks with a 
Mach number M < 1.3, Ryu et al.2003); however in 
the following we neglect this third condition, in order 
to have a better comparison with the results obtained 
with the VJ method (see below). 

The shock discontinuity in the simulation is typi- 
cally spread over a few cells, thus following Hallman 
et al.(2004) for each patch of candidate shocked cells 
we define the shock center with the position of the cell 
in the shocked region where V • v is minimum and 
calculate the Mach number of the shock from Eq]4j 
where T2 and T\ are the post and pre-shock tempera- 
ture across the shock region. The Mach numbers mea- 
sured along the three coordinate axes are finally com- 
bined to compute the Mach number of the shocked cell: 
M = (Ml + Ml + Ml) 1 / 2 . 

More specifically, in the case the shock-jump is as- 
sumed to happen in 1 cell, T2 is the temperature of the 



shock center, while in the case that the jump is spread 
over 2, 4, 6, . . . 2 n cells T2 and T\ are the tempera- 
tures of the two cells at distance of n cells (in opposite 
direction) from the center of the shock. 

In the following we will refer to this method as the 
77 method. 



5.3 The Velocity Jump Method. 

Conservation of momentum in the reference frame of 
the shock yields: 



P\v\ = P2V2, 



(6) 



with the same notation used in EqsfJHU In the ideal 
case in which the pre-shocked medium is at rest and 
in thermal and pressure equilibrium, the passage of a 
shock with velocity v s leaves a clear imprint as a ve- 
locity difference, Av, between the shocked and pre- 
shocked cells. The relationship between Av and Mach 
number in the case of hydrodynamical shocks can be 
obtained by combining Eqj6] with Eq[3] and by trans- 
forming the velocities from the shock frame to the Lab 
frame : 



3 1 



M 2 



(7) 



where v s = Mc s and c s is the sound velocity com- 
puted in the pre-shocked cell. 

The procedure we adopt to identify shocks and re- 
construct their Mach numbers is the following : 

• we consider candidate shocked cells those with 
V • v < (calculated as 3-dimensional velocity diver- 
gence to avoid confusion with spurious 1-dimensional 
compressions that may happen in very rarefied envi- 
ronments); 

• since shocks in the simulation are typically spread 
over a few cells, as in the case of the TJ method, we 
define the shock center with the position of the cell in 
the shocked region with the minimum divergence; 

• we scan the three Cartesian axes with a one- 
dimensional procedure measuring the velocity jump, 
^■Vx,y,z, between a few cells across the shock center. 
In the case the shock-jump is assumed to happen in 1 
cell Av xy ,z is calculated between the shock center and 
the pre-shocked cell, while in the case that the jump is 
spread over 2, 4, 6, . . . 2 n cells Av X:VjZ is calculated 
between two cells at distance of n cells (in opposite di- 
rection) from the center of the shock; 

• the Mach number of the shock is given by EqJTJ 
where the sound speed is that of the pre-shock region 
(the cell with the minimum temperature); 

• we finally assign to shocked cells a Mach number 
M = (Ml + M 2 + Ml) 1/2 , that minimizes projection 
effects in the case of diagonal shocks, and restrict to 
shocks with M > 1. 

In the following we refer to this procedure as the ve- 
locity jump (VJ) method. 
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Figure 8. Scatter plots for the value of Mach numbers from Monte Carlo extracted populations of non shocked cells, in filaments, outskirts and clusters (see 
text). Temperature fluctuations are extracted from the corresponding distributions in Fig[7] The red curves give the exact solution from Eq|4] 





Figure 9. Scatter plots for the value of Mach numbers from Monte Carlo extracted populations of non shocked cells, in filaments, outskirts and clusters (see 
text.) Velocity fluctuations are extracted from the corresponding distributions in Fig|7] The red curves give the solution from Eq|7] 



6 UNCERTAINTIES IN SHOCK DETECTING SCHEMES 

In this Section we discuss the uncertainties of the meth- 
ods presented in the previous Section, discuss the effect 
of the re-ionization on the characterization of cosmo- 
logical shocks and briefly compare results from the VJ 
and TJ approaches. 



6.1 Reconstruction of the shock discontinuity 

Although the shock discontinuity in ENZO is found 
to be well reconstructed within 2-4 cells (e.g., Tasker 
et al. 2008), the risk that comes from the application 
of procedures based on cell-to-cell jumps (or jumps 
between a few cells) is to underestimate the value of 
the Mach number of the shock. We performed several 
shock-tube tests with ENZO with the same numerical 
setup used in the cosmological simulations, in order 
to evaluate the convergence of the value of the shock 
Mach number with the number of cells used to cal- 
culate jumps. We find that a reasonable convergence, 
within 10-30% for M < 10, is already obtained with 
the VJ method in the case that the velocity jump is eval- 
uated across three cells (n = 1, where n is the distance, 



in term of cells, between the shock center and the pre or 
post-shock cells, Sects. 5.2-5.3), and that convergence 
is reached for n ^ 2. On the other hand, the velocity 
pattern in cosmological simulations is complex and the 
risk of procedures based on jumps evaluated with large 
n in our simulations is to mix together signals produced 
by different shocks, and also to be affected by gradients 
in thermodynamic al variables due to clumps of bary- 
onic matter produced by the process of structure for- 
mation. In Fig,[6]we show a map of the Mach numbers 
obtained with the VJ method for a galaxy cluster in the 
AD125 run, by assuming a cell-to-cell (two cells) ve- 
locity jump, and n = 1, n = 2 and n = 4 jumps. It is 
clear that for n ^ 2 (jumps based on ^5 cells, ^625 
kpc) different shock-patterns and clumps of gas matter 
start to be mixed together and shocks become poorly 
characterized. 



Similar results are found in the case of the TJ 
method, thus we conclude that reconstructing the shock 
discontinuities in our numerical simulations with n = 
1 (jumps based on 3 cells) provides the best compro- 
mise. 
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Figure 10. Maps of a 20Mpc region centered on a Mtot ~ 1O 14 M0 galaxy cluster, with 125fcpc of cell resolution (the line of sight width is one cell). Left 
column shows the baryon density (top panel) and baryon temperature (bottom panel), in this case taken from the same region of the C0125 run); the central 
and the right columns show the maps of Mach number according to the VJ scheme (bottom panels) and to the TJ scheme (bottom panels). Shocks are shown 
both in the case of no-reionization (AD125, central column) and for the case of re-ionization (C0125, right column). 



6.2 Uncertainties in the TJ and VJ methods 

As already pointed out in Sectj5J a major limitation of 
the analysis of shocks in post processing comes from 
the fact that the dynamics and thermodynamics of the 
gas in the simulations is more complex than in the ideal 
case in which Eqsf3]-f5]and|7]are derived. In this sub- 
section we discuss the uncertainties that come out with- 
out including the modeling of re-ionization in our pro- 
cedure. 



6.2.1 TJ method 

The temperature distribution in simulations is very 
complex and temperature gradients across non- 
shocked cells are common by-products of the struc- 
ture formation process. The passage of a shock in a 
medium with a complex temperature distribution par- 
tially smooths out pre-existing gradients in the thermo- 



dynamical variables and creates new shock-induced 
discontinuities. 

One possibility to estimate the level of uncertain- 
ties in the application of the TJ method is to introduce a 
passive modification of the post-shock temperature in 
Eqj4]according to the value of a typical cell to cell tem- 
perature jumps across non shocked cells, and compare 
the resulting Mach number with that from EqQ] in its 
original form. Obviously this procedure assumes that 
these jumps are representative of pre-existing temper- 
ature gradients, where shocks are presently found, still 
there is no clear argument for which this unavoidable 
assumption is not statistically reasonable. 

We consider as non-shocked cells those that do not 
satisfy, at the same time, the TJ and VJ criteria for 
shocked cells, and extract the values of their cell to cell 
temperature jumps, ST, in different cosmic environ- 
ments from the AD 125 simulation at z = 0. To follow a 
very conservative procedure we consider only temper- 
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ature jumps across a sub-sample of non shocked cells 
that are at least three cells far away from any shocked 
cell. 

We characterize the cosmic environment by means 
of the total matter density in cells : 

• 0.01 ^ Ptot/Pcr < 3: voids and under-dense re- 
gions, 

• 3 ^ ptot/pcr < 30 : filaments and sheets, 

• 30 ^ Ptot/Pcr < 50: cluster outskirts, 

• Ptot/Pcr ^ 50: galaxy clusters. 

where p tot = p + pdm is the total matter density 
and p cr is the critical density of the universe. These are 
expected to mark the different kind of structures of the 
cosmic web (e.g. Dolag et al.2006; Shen et al.2006). 
In Fig|7] we plot the differential distribution of (the 
module of) temperature jumps across the considered 
sub-sample of non-shocked cells, ST, normalized to 
a reference value of the local temperature, \fT% ■ T\, 
for the different density regimes. The peak of the 
distribution is found at ST/^T 2 ■ T\ sa 0.5. In the 
case of filaments and voids the distributions extend 
at larger values, although large temperature scatters, 
ST/ \JT~2 • T\ « 10, are only found in voids, where they 
represent less then 1 percent of all cellfl 

extremely rare in the case of filaments and are 
found for only a few percent of the cells in the voids. 
For the values of T\ representative of clusters, outskirts 
and filaments in our simulation, we allow T 2 to vary 
and run Monte Carlo extractions of ST extracted across 
non-shocked cells with temperature T 2 in the same en- 
vironment. We then compared the estimate of the shock 
Mach number by Eqf4]with the one obtained by using: 



T 2 ±\ST\ (5M 2 - 1)(M 2 + 3) 



Ti 



16M 2 



(8) 



Figure 7 shows the typical scatter introduced in 
the T 2 /T\ vs M plane by the presence of realistic 
(i.e. measured in non-shocked cells in our simulations) 
pre-shock fluctuations in the temperature, for differ- 
ent cosmic environments. The red line shows the ideal 
case: given a ratio T 2 /T\ the degree of uncertainty on 
M due to the presence of pre-shock fluctuations in the 
simulation can be grossly evaluated by an horizontal 
cut across the distribution of the data points. This scat- 
ter increases as the Mach number decreases and, at a 
given Mach number, it is typically smaller in environ- 
ments with larger over-density. 



6.2.2 V J method 

Complex velocity fields arise naturally during the for- 
mation of virialised structures in simulations (Bryan & 
Norman 1998, Sunyaev; Norman & Bryan 2003; Dolag 
et al. 2005; Vazza et al. 2006; Iapichino & Niemeyer 

3 We also note that rare large temperature jumps in low density environ- 
ments, not associated with shocks, might come from numerical artifacts in 
regions where gas flows are supersonic and the computation of the internal 
energy is more subject to numerical uncertainties. 



no reionization 
H&M reionization 
postproc T= 1 00 
postproc T=100Q 
postproc T= 1 A 4 




10000 



Figure 11. Distribution of kinetic energy flux in shocks according to the 
VJ method, for a cubic volume of side AOMpc and resolution 125kpc. 
Curves are drawn for the case without re-ionization (black solid), for the 
Haardt & Madau (1999) re-ionization scheme applied in post-processing 
(blue dashed) and for different choices of a fixed Tfi oor temperature floor 
(color coding is labeled in the panel). 



2008) that however are expected to be smaller than the 
velocity jumps driven by the passage of a shock across 
the same regions. A more complex situation is that 
of non virialised structures where laminar flows may 
produce relatively strong velocity gradients across the 
cells. An example is given in Fig. where we report 
the differential distributions of the velocity gradients, 
Sv, normalized to the maximum value of the sound 
speed in each pair of cells, obtained for the same sub- 
sample of non-shocked cells considered in the previous 
sub-Section. The distributions peak at Sv/c s rj 0.5, al- 
though tails extending towards larger values are found 
in the distributions of voids and filaments. These tails 
are mostly due to velocity gradients measured across 
accelerated laminar flows (where the kinetic energy of 
the gas may become even larger than the thermal en- 
ergy) that move from low to higher density regions and 
are present in a small fraction of the volume of fila- 
ments and voids. 

In order to grossly estimate the strength of the uncer- 
tainties in the case of clusters, outskirts and fila ments, 
we follow a procedure similar to that in Sect J6.2.ll 
For these different environments we fixed a value of 
Av/c s , run a Monte Carlo extraction of Svjc s from 
non shocked cells in the simulations (FigjT) and for 
each trial calculated the Mach number from : 



A , . 3 l-M 2 

Av ± \Sv\ = 

1 1 4 M 



1 ± 



Sc, 



(9) 



This equation accounts for both pre-shock gradi- 
ents in the velocity and in the sound speed across non 
shocked cells. Gradients in c s are driven by gradients 
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Figure 12. Distributions of shock Mach numbers at different cosmic environments, for a (80Afpc) 3 volume of the ADl2b run with post-processing re- 
ionization (left) and for the C0125 run (right), with the VJ method (solid lines) and with the TJ method (dot-dashed lines). 



in the temperature distribution of the cells and are eval- 
uated by a Montecarlo extraction of the temperature 
variations in Fig. [7] 

In Figj9]we report Av/c s vs M from our Monte 
Carlo extraction compared to the calculations in the 
ideal case (EqJT]). This result should be compared with 
that in Fig JS] obtained for the TJ method and the degree 
of uncertainty on M can be grossly evaluated by an 
horizontal cut across the distribution of the data points. 
As expected, in the case of clusters and outskirts the 
scatter in the two cases is quite similar, although in the 
case of weak-moderate shocks crossing filaments and 
outskirts the scatter is less pronounced than that of the 
TJ method (Fig©. 

6.3 Modeling the re-ionization. 

The role of re-ionization is of primary importance to 
study the properties of shocks outside galaxy clusters. 
In adiabatic simulations, regions far away from intense 
structure formation are very cold due to the lack of in- 
tense shock heating and due to cosmic expansion. In 
these environments, any additional source of heating 
(such as re-ionizing radiation from AGN and/or mas- 
sive stars feedback) causes a dramatic increase of the 
local temperature and sound speed. Thus the tempera- 
ture distribution across cells in these regions is strongly 
affected by the modeling of the re-ionization in the 
simulations, and this implies an additional uncertainty 
in the characterization of shocks. This is expected to 
be particularly relevant in all shock detecting schemes 
where temperature plays a role. 

Therefore in this Section we highlight the main ef- 
fect of cosmic re-ionization on shocks FigfTO] shows 
the maps of the detected shocks in a 20Mpe cubic re- 
gion extracted from the AD125 simulation and cen- 
tered on a M to t ~ 10 14 M Q cluster. Results are re- 
ported, by calculating shock-jumps across three cells 
(n=l, Sect. 5.2, 5.3), in the case of no re-ionization 



and of a Haardt & Madau (1996) re-ionization scheme. 
As expected the Mach number of shocks decreases in 
simulations with re-ionization due to the increase of 
the sound speed produced by the re-ionization back- 
ground. This effect is stronger in the cold outermost 
regions, while the properties of cosmological shocks 
in galaxy clusters are not affected by the re-ionization 
background. Re-ionization also allows to better de- 
scribe shocks around filamentary structures in low den- 
sity regions that are not seen in the case of simulations 
without re-ionization (FigfTOl). This is because without 
re-ionization these regions have temperature so small 
that the temperature floor (IK) adopted (in the outputs) 
by ENZO artificially affects their temperature distribu- 
tion. 



In FigfJTJwe report the kinetic energy flux through 
shocked cells as calculated by means of the VJ method. 
The kinetic energy flux, Eki n = is reported for 

different numerical modeling of the re-ionization: three 
different temperature floors (W^K, 10 3 and 10 2 K), 
Haardt & Madau (1999) model, and no re-ionization. 
We find that a fixed temperature floors, which is cus- 
tomary used in several papers to mimic the effect of 
re-ionization (e.g. Ryu et al. 2003), produces some ar- 
tificial piling up or flattening in the distributions of the 
energy flux through shocks at large Mach numbers. 
This is because the temperature background, Tfi oor , 
changes artificially the speed of the sound in environ- 
ments with lower temperature and the Mach number of 
shocks in these environments is affected by Tp oor , de- 
creasing artificially with increasing Tfi oor . This further 
supports the requirement of a physically meaningful 
treatment of re-ionization in a post processing proce- 
dure. As already discussed, the post processing fitting 
procedure described in Sect[4]closely resembles the ef- 
fect of the physically based Haardt & Madau (1999) 
re-ionization scheme. 
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Figure 13. Maps of a region of side 80Mpc from the AD125 run with post-processing re-ionization, showing gas temperature (left) and Mach number 
measured with the VJ method (right). The width along the line of sight is 125fcpc. 



6.4 Basic Comparison between VJ and TJ methods 

In this Section we briefly compare the results obtained 
from the VJ and TJ approaches. Here we focus on re- 
sults obtained with our fiducial numerical treatment of 
the re-ionization as this treatment will be used in the 
following of the paper; a more detailed comparison is 
ongoing and will be presented in a forthcoming paper 
(Vazza et al. in prep). 

In the ideal case the two approaches must select the 
same population of shocked cells. In reality we find 
that, when shock-jumps are calculated across 3 cells, 
about 85 per cent of the shocked cells in our simula- 
tions are selected at the same time by the conditions in 
the VJ and TJ approaches. In the case of clusters and 
cluster outskirts the two approaches select the same 
population of shocked cells, while these differences 
typically arise from shocks in low temperature regions. 

In the case of clusters and outskirts the velocity 
and temperature variations across non shocked cells are 
relatively small (Fig. [7) and this allows constraining the 
Mach number of shocks by means of both the TJ and 
VJ approach. Still the statistical uncertainties for weak 
shocks with the TJ method are expected to be slightly 
larger than those with the VJ (Figsj8]&[9]). 

A comparison between the statistical description of 
the properties of the shocks with VJ and TJ approaches 
is shown in FigfT2] where we report the Mach number 
distribution of shocks extracted from the AD 125 run 
with our post-processing treatment of re-ionization and 
from the corresponding CO 125 with run-time treat- 
ment of re-ionization. Basically the results from the 
two methods are fairly similar in the case of clusters 
and cluster outskirts, and no remarkable differences are 



found also in the case of filaments and voids. This sug- 
gests the important point that the characterization of 
shocks in these environments is statistically solid, as 
two independent approaches lead to basically similar 
results. We also note that the Mach number distribu- 
tions extracted from simulations with post processing 
re-ionization are similar to those produced with run 
time re-ionization, and this further suggests that our 
treatment of re-ionization in post processing allows sta- 
tistically valid studies of shocks. 

In the next Sections we shall use the VJ method 
to study shocks properties. This is because we believe 
that in the case of weak-moderate shocks, especially in 
lower density regions, the VJ is less affected by uncer- 
tainties (e.g., Figs. 8-9). 



7 RESULTS 

In this Section we present the main results obtained for 
the full set of simulations by making use of the VJ 
method and by calculating shock-jumps across three 
cells, i.e. n = 1 (unless specified). 



7.1 Detected shocks and Maps. 

Shocks fill the simulated volume in a very complex 
way (e.g. Miniati et al.2001, Ryu et al.2003). In FigfT3] 
we show a 125 kpc cut of a cubic region of side 80Mpc 
from the AD125 run at z = with post-processing 
re-ionization, reporting gas temperature and detected 
shocks with Mach numbers reported in color code. 

We find that 10 — 20 per cent of the cells in the 
simulated volume host shocks at present epoch, with 
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Figure 14. Distribution of shocks Mach number for the whole simulated 
volume of the AD125 with post-processing re-ionization, for different cos- 
mic environments. Dot-dashed lines show the distributions obtained with 
velocity jumps evaluated across three cells (n = 1), while dashed lines 
shows distributions obtained with cell-to-cell velocity jumps (ra = 0). 



Time Evolution 
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Figure 15. Time evolution of the Mach number distribution for the C0125 
run, from z = 7.2 to z = 0.0. Only a sub sample of redshifts is shown for 
clarity. 

the percentage of shocked cells increasing in denser en- 
vironments. Filamentary and sheet-like shocks pattern 
are usually hosted in low density regions and at the in- 
terface of filaments, following the shape of the cosmic 
web. This kind of shocks follows the first infall of bary- 
onic matter onto accreting structures, and generates an 



abrupt increase in temperature due to the jump from 
a re-ionization dominated temperature to the gravita- 
tionally dominated one. These shocks are commonly 
defined as "external shocks" (Miniati et al.2001), and 
they are the strongest within the simulation, having 
M » 10. Shocks surrounding galaxy clusters form 
spherically shaped boundaries at a typical distance of 
about 2R v i r from the cluster center, while shocks mov- 
ing inwards the virializing region are found more irreg- 
ular and weak, with M < 3. This kind of shock waves 
is commonly defined as "internal shocks " (Miniati et 
al.2001). Slightly stronger shocks (i.e. M ~ 3) inside 
R v i r are episodically found in our simulations in case 
of merger events. In this case the violent relaxation 
due to the fluctuation of the gravitational potential may 
cause infall of the pre-shocked gas onto the shock dis- 
continuity increasing the Mach number (Springel & 
Farrar 2007), while other kind of strong shocks are re- 
verse shocks propagating trough the innermost regions 
of accreting and cold sub clumps, which keep them- 
selves at the pre-shock virial temperature for several 
Gyrs during their orbiting around the center of the main 
cluster (Tormen, Moscardini & Yoshida 2004). 

An issue which is poorly addressed in the literature 
is the distribution function of shocks with their Mach 
number. Fig{l4] shows the Mach number distribution 
of the shocks detected in our total l45Mpc cubic vol- 
ume at present cosmic epoch. This figure also shows 
the effect of using cell-to-cell velocity jumps to recon- 
struct the Mach number of shocks with the VJ scheme. 
The number of stronger shocks, that are well recon- 
structed within 3-4 cells, increases with n = 1 and this 
produces a flattening of the differential distribution of 
shocks inside the simulated volume. As shown in Sec- 
tion 6.1, using a larger number of cells to reconstruct 
the Mach number does not improve the characteriza- 
tion of shocks, yet the risk becomes to mix different 
shock patterns and sub-structures in the simulations. 

The overall differential distribution of shocks with 
their Mach number in the cosmic volume is very steep, 
with a ~ -1.6 (with MdN/dM oc M a ), and the bulk 
of the detected shocks at any Mach number is found in 
the low density regions, which fill the majority of the 
volume in the simulations. The Mach number distribu- 
tion of detected shocks becomes increasingly steeper 
moving towards dense environments: a « —3 to —4 is 
found in clusters and their outskirts. 

The evolution with time of the differential Mach 
number distribution is given in FigJT5] for the case of 
a 80 Mpc cubic region extracted from the C0125 sim- 
ulation, which was designed to have a suitable time- 
sampling in the analysis of outputs. We find that be- 
fore the epoch of re-ionization, z > 6, roughly 30% 
of the simulated volume is shocked. Then as soon as 
re-ionization plays a role, the temperature of the simu- 
lated volume increases and the Mach number distribu- 
tion of shocks at redshift z ~ 3 — 6 undergoes a dra- 
matic change becoming very steep and dominated by 
weak shocks. With decreasing redshifts, temperature in 
low density regions gradually decreases and the Mach 
number distribution becomes gradually flatter with the 
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Figure 16. Left:3-D rendering of baryon density for a cubic region of side 80Mpc, for the C0125 re-simulation at z = 0. Color coding goes from dark 
blue (p ~ 10 — 31 gr /cm 3 ) to pale blue (p > 10 — 29 gr/cm 3 ). Right: thermalized energy flux through shocks, for a slice of depth V2bkpc and centered to 
encompass the two massive merging clusters shown in the left panel (A and B). 



fraction of shocked cells reaching ~ 15 per cent at 
present epoch. 



7.2 Energy Flux and thermalised energy 

The energy flux converted into thermal energy of the 
gas at a shock is given by the Rankine-Hugoniot jumps 
conditions, which relate the flux of the kinetic energy 
crossing the shock, E^in, and the resulting thermal flux 
in the post-shock region, fa. This relation can be ex- 
pressed by means of a simple 5(M) parameter (e.g. 
Ryu et al.2003): 



5{M) = fa/ h = v 2 



Eth,2 



in,l 



(10) 



where E t ^\ and E t h,2 are the thermal energies in 
the pre- and post-shock regions, E^inx is the kinetic 
energy of the shock, and T is the adiabatic exponent 
(T = 5/3). It is useful to express 5(M) by means of 
the Mach number (e.g. Kang et al.2007): 



S(M) 



r(r - l)M 2 R 



2TM 2 - T + 1 
f+1 



(11) 



where R is the density compression factor: 



R = ^ 
Pi 



r + i 



r - i + 2/M 2 



(12) 



We notice that Eq[l0] strictly holds only in case of 
a negligible CR energy density, otherwise the feedback 
of these CR on the shock itself is expected to severely 
decrease the efficiency of thermalisation of the kinetic 
energy flux (see next Section). 

FigHH (Right panel) shows a 2-dimensional cut, 
with depths 125 kpc, of the measured thermal energy 
flux in shocked cells, at the present epoch and for a re- 
gion centered around two massive (M ~ 4 • 1O 14 M 
and M ~ 1O 15 M ) galaxy clusters. These clusters be- 
long to a large scale filament (see Left panel of FigJT6l>, 
for which we provide also a 3-dimensional rendering 
of the thermal energy flux through shocks (Left panel 
ofFiginp 

The differential distribution of the thermalised flux 
as a function of the Mach number of the shocks is 
reported in Fig{T8] This shows the differential distri- 
bution calculated for the total volume at the present 
epoch and normalized to the comoving volume of 



4 We generated 3-dimensional distribution of data by means of the visual- 
ization tool VISIVO (Comparato et al.2007, http://visivo.cineca.it I 
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Figure 17. Left:3-D rendering of the dissipated energy flux for the same region as in Fig l 161 Color coding goes from blue (f t h ~ 10 33 erg/s) to yellow 
(fth ~ 10 38 erg/s) to red {fth > 10 41 ercj/s). Right: energy ratio between injected CR energy flux and thermal energy flux in shock waves, for the same 
slice of right panel of Fig l 161 



(IMpc/h) 3 . Solid lines give average values, while 
dashed shadows give the variance spanned by the six 
different 80Mpc cubic sub samples of the AD 125 sim- 
ulation. The variance in the distribution is fairly small, 
~ 30 per cent at the peak, although this increases for 
the stronger shocks that are rare. 
We find that the total processed thermal energy across 
cosmological shocks in our simulations is f t ^ « 4 • 
10 47 ergs/ s at the present epoch. This is of the same 
order of magnitude of the value of the total processed 
thermal energy found by Ryu et al.(2003) and by 
Pfrommer et al.(2006), for the same ^ 145 M pc cubic 
volume. However, as discussed in Sect j3.2| the deficit 
in massive halos in our clusters sample may induce a 
slightly smaller level of thermalised energy flux in the 
volume. 

For Mach numbers < 20 (which provide about the 99 
percent of the total thermal flux in the simulated vol- 
ume) the distribution in Fig{l8]has a tri ~ —2.7 (with 
fth(M)M oc M ath ), and is steeper than that in Ryu et 
al.(2003), a t h ~ —2, while is consistent with that in 
Pfrommer et al. (2006), a th « -2.5. 
We find that « 70 per cent of the total thermal en- 
ergy flux dissipated at shocks comes from the virial 
region of galaxy clusters (because of their large mat- 
ter density) and that the bulk of the thermalisation hap- 



pens at shocks with M « 2 (Fig [18]). These relatively 
weak shocks are also responsible for the bulk of the 
thermalisation in lower density environments, although 
stronger shocks may provide a sizable contribution in 
these regions. 



The time evolution of the distribution of the ther- 
mal energy dissipated at shocks as a function of the 
shock-Mach number is an important issue. A rele- 
vant example is reported in Fig{l9] that is obtained 
for the same volume of Fig. [15] (C0125). The evo- 
lution of the overall distribution follows a similar be- 
havior with cosmic time found for the number distri- 
bution of shocks, with strong shocks becoming more 
frequent at evolved times when the energy density of 
the background becomes lower (see also Pfrommer et 
al. 2006). The integrated (over cosmic time) thermal 
energy dissipated at shocks in our (145Mpc) 3 volume 
is Eth ~ 2 • W 64 ergs, which is consistent with the 
values reported in Pfrommer et al. (2006) and Ryu et 
al.(2003), also by taking into accou nt th e deficit in the 
halos mass in our simulations (Sect l3.2b . 



16 F. Vazza, G. Brunetti, C. Gheller 



AT 2 5 + i"ei o r i 7 a t " o i"i 



voids 




ID 'S§ 1D0G 



Figure 18. Distribution of the thermalized energy flux at different over- 
densities, for the whole AD125 and normalized to a comoving volume of 
(IMpc/h) 3 . The shadowed regions show the cosmic variance within our 
6 simulations, the dot-dashed line shows the global average within the sam- 
ple. 



7.3 Acceleration of Cosmic Rays 

The injection and acceleration of Cosmic Rays at 
shocks is a complex process. It is customary to describe 
the acceleration according to the diffusive shock accel- 
eration (DSA) theory (e.g. Drury & Voelk 1981; Bland- 
ford & Ostriker 1978). This theory applies when parti- 
cles can be described by a diffusion-convection equa- 
tion across the shock. There is some general agreement 
on the fact that strong shocks may channel a substantial 
fraction of their energy flux into the acceleration of CR 
which in turn should back react modifying the structure 
of shocks themselves. Recent advances rely on the the- 
ory of non linear shock acceleration, which describes 
the acceleration of CR in shocks whose structure is 
modified by the back-reaction of CR energy (e.g., Elli- 
son, Baring, Jones 1995; Malkov 1997; Kang, Jones & 
Gieseler 2002; Blasi 2002, 2004a; Kang & Jones 2005; 
Amato & Blasi 2006). The most relevant uncertainty 
in the description of the particle acceleration at these 
shocks is the injection model, i.e. the probability that 
supra-thermal particles at a given velocity can leak up- 
stream across the sub shock and get injected in the CR 
population. This is because even a small variation of 
the injection momentum, Pi n j, of supra-thermal parti- 
cles produces a large difference in the estimate of the 
injection efficiency at shocks (e.g. Blasi 2004b). An 
other major hidden ingredient is the amplification of 
the magnetic field (perpendicular component) down- 
stream due to CR driven instabilities and adiabatic 
compression, as this magnetic field self-regulates the 
diffusion process of CR and supra-thermal particles 
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Figure 19. Time evolution of the distribution of thermal energy flux at 
shocks for the same volume as in Fig[l5] from z = 6.1 to z = 0.0. Only a 
sub sample of redshifts is shown for clarity. 

(i.e. the Larmor radius) and considerably affects the 
value of pi n j. 

An additional difficulty which comes out is that a post- 
processing approach, as that followed in our paper, 
does not allow us to account for the dynamical con- 
tribution of CR accelerated at cosmological shockfl 

With all these caveats in mind, we follow the 
approach adopted by Ryu et al.(2003) in which the 
thermalisation is calculated by means of the standard 
Eqs flQl - flTl and the CR acceleration at shocks is calcu- 
lated by making use of numerical results of non linear 
shock acceleration which adopt a numerical descrip- 
tion of the thermal leakage to model the injection of 
particles in the population of CR upstream (Kang & 
Jones 2002, KJ02). These numerical results provide an 
estimate of the ratio between the energy flux trough a 
shock and the energy flux which is channeled into CR 
acceleration at the shock by means of a simple param- 
eter, r](M) = fen/ fd>, which depends on the Mach 
number of that shock. 

FigfTT] (Right panel) maps the ratio between CR 
and thermal energy flux for the same region reported 
in Right panel of FigJT6j This clearly shows the role 
played by the Mach number in setting the level of 
the injection of CR in the various environments. Since 
the ratio r](M)/5(M) (= fen/fth) increases with 
the Mach number of the shocks, the highest values of 
Icr/ fth ^ found in low density regions, at the in- 
terface layers of filaments or in the outermost regions 
of galaxy clusters, where a substantial population of 
relatively strong shocks is present. On the other hand 

5 Attempts to model this dynamical contribution in cosmological simula- 
tions have been recently developed (Pfrommer et al. 2006) 
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the lower values are typically found in galaxy clusters, 
where the Mach number distribution is steep and strong 
shocks are rare. 

The distribution of the energy flux injected in CR 
as a function of Mach number is reported in Figj20l(Zeff 
panel); this refers to the total simulated volume at the 
present epoch. We find that the bulk of the CR accel- 
eration at present epoch takes place in galaxy clusters, 
however also filaments are expected to contribute sig- 
nificantly to the overall acceleration process. The over- 
all distribution has a well defined peak which is an- 
chored at M 2 and has a slope (between M ~ 
2 - 20) a C R ps -2 (with f CR {M)M oc M acR ). 

The value of the Mach number at the peak is 
close to (even if slightly smaller than) that found by 
Ryu et al.(2003) (M ~ 3), while the distribution is 
steeper than that reported in Ryu et al.(2003) (where 
a CR ~ —1-5). Since we use an approach equivalent 
to that in Ryu et al.(2003) to evaluate the CR acceler- 
ation, this difference is likely related to our different 
shock detecting scheme, to our different modeling of 
the re-ionization process and also to the slightly dif- 
ferent value of the eg parameter (erg = 0.8 in Ryu et 
al.2003; we discuss in the Appendix the influence of erg 
in the shock statistics). A comparison with the results 
in Pfrommer et al.(2006) is more difficult since these 
authors use a Lagrangian Smoothed Particles Hydro- 
dynamics code which also include CR dynamics and 
a completely different approach in the calculation of 
the CR injection at shocks. The overall distribution of 
the energy flux injected in CR reported in Pfrommer et 
al.(2006) has a slope acR ~ —1-8 and is actually in 
between our results and those of Ryu et al. (2003). 

For seek of completeness, in Figf20] (Right panel) 
we also report the overall distribution of the energy flux 
injected in CR by adopting the injection efficiency of 
CR at modified shocks by Kang & Jones 2007 (KJ07). 
These recent calculations account for the Alfven wave 
drift and dissipation in the shock precursor and this 
yields a value of rj(M) which is smaller than that 
adopted by Ryu et al.(2003) (at least for M < 20). 
As a consequence the resulting distribution of the en- 
ergy flux dissipated at shocks into CR acceleration as a 
function of the shock-Mach number is flatter than that 
obtained by adopting KJ02 and « 50 per cent less en- 
ergy is expected to be channeled into CR acceleration. 

FigfJT] shows the evolution with time of the ratio 
fcR/fth f° r the same volume (C0125) considered in 
Figs. 1151 and [T9l The value of fc R /fth as measured 
at the present epoch, fcR/fth ~ 0-2, is a factor «2 
smaller than that found in Ryu et al.(2003). By adopt- 
ing the injection efficiency of CR of KJ07 the ratio 
fcR/fth is even smaller, about 0.1, as the injection ef- 
ficiency of CR in KJ07 is smaller than in KJ02 (at least 
for shocks with M < 20). Although Fig.HTIshows that 
CR in dense regions do not provide a relevant back- 
reaction on the thermal pool during their acceleration 
(this justifies the use of Eqs. 10-12 in these environ- 
ments), the larger values of /cr/ fth in low density re- 
gions and at early cosmic times suggest that follow- 
ing (run-time) the dynamics of CR and the (self con- 
sistent) non-linear shock thermalisation and CR accel- 
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Figure 21. Evolution with time of the total injection ratio fen/ fxH for the 
same sub sample as in Fig. 15 and !19l shown for different environments. The 
upper panel is for the acceleration model outlined in Kang & Jones (2002), 
while the lower panel is for the acceleration model of Kang & Jones (2007). 

eration is mandatory in future studies with Eulerian- 
cosmological simulations. 

7.4 Shocks in Galaxy Clusters. 

In this Section we focus on the shock statistics and CR 
injection in galaxy clusters and briefly discuss their de- 
pendence on the cluster dynamics. We study shocks in 
four representative massive galaxy clusters extracted 
from the AD125 simulation, at z = 0: 

• CI: aMtat ~ 7-lO 14 M cluster in a relaxed state; 

• C2: a Mtat ~ 7 • 1O 14 M cluster subject to an 
ongoing minor merger with a sub clump with mass 
Mtot ~ 2 • 1O 13 M ; 

• C3: a M to t ~ 1 • 1O 15 M cluster approaching a 
major merger with zero impact parameter, with a com- 
panion cluster (with M to t ~ 4 • 1O 14 M ) that is at a 
distance of ~ 1.3R v i r from the main cluster center; 

• C4: a M^t ~ 7.5 • 1O 14 M cluster in a post- 
merging phase (the merger occurs in the simulation ~ 
2 Gyr in look back time). 

These clusters are reported in Figf22l that shows 
maps of projected baryon density in a (4i?„j r ) 2 regions 
centered on clusters, and maps of the Mach number 
measured in slices crossing the same regions. 
In the case of C2 (minor merger) and C3 (major 
merger) relatively weak, M « 2 — 2.5, shocks are 
found inside R v i r , while in the case of C4 (post- 
merger) merger shocks have already moved outside the 
internal region of the cluster, and their strength is in- 
creased as the ambient temperature in cluster outskirts 
decreases. 
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Figure 20. Distribution of the injected CR flux at different over-densities, for the whole AD125 run with post-processing re-ionization. The shadowed regions 
show the cosmic variance within our 6 simulations. Left panel shows the measured distribution according to the KJ02 recipe for the CR injection, while right 
panel is for the KJ07 recipe. 




Figure 22. Left panels: maps of projected baryon density for the 4 galaxy clusters of Sect l7.4l Every map has a depth along the line of sight of twice the virial 
radius of the correspondent cluster.Ri'g/if panels: slabs of 125 kpc along the line of sight, showing the maps of Mach number for the same objects as in left 
panels. For displaying purposes, pixels in the images have been oversampled and convolved with a Gaussian kernel with a FWHM=2 cells. 
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Figure 23. Volume averaged profile of the mean Mach number of shocks 
for the galaxy clusters in Fig |22l 



8 DISCUSSION AND CONCLUSION. 

In this paper we have investigated the properties of 
shocks in cosmological numerical simulations. This 
subject is particularly intriguing as shock waves prop- 
agating trough LSS are the responsible for the heat- 
ing of the ICM and may be important sources of CR 
in the Universe. This subject has been already inves- 
tigated in several papers under different numerical ap- 
proaches (Miniati et al.2000; Miniati et al.2001 ; Keshet 
et al.2003; Ryu et al.2003; Pfrommer et al.2006; Kang 
et al.2007; Pfrommer et al.200'®. 

We study shock waves by means of a post pro- 
cessing procedure. Although this is similar to Ryu et 
al.(2003) and Kang et al.(2007), our approach differs 
from previous ones in several points. First we use a 
different numerical codes, the public version of ENZO 
(e.g. Bryan & Norman 1997), to simulate LSS (Sec©. 
Second we adopt a more appropriate treatment of the 
re-ionization in our simulations (Secj4j) and third we 
use a different approach to catch shocks in our sim- 
ulations and to measure their strength (VJ method, 
SectEl. 



The volume averaged Mach number of shocks in 
the four galaxy clusters is reported as a function of dis- 
tance from cluster centers in Fig[23] Within the virial 
radius shocks are weak in line with expectations from 
semi-analytical models that indeed found shocks with 
M > 3 extremely rare in galaxy clusters (Gabici & 
Blasi 2003). This is also highlighted in Fig[24l where 
we show the distribution of the thermal flux dissipated 
at shocks with shock-Mach number; distributions in 
different clusters are reported normalized to the vol- 
ume of the most massive cluster. All distributions are 
steep, with differences from cluster to cluster due to 
the effect of the dynamical state of the clusters. In- 
side R v i r CI, C2 and C3 have similar distributions, 
while C4 shows some excess of rare shocked cells with 
M f« 3 — 7. On the other hand, an excess of shocked 
cells with M 3 — 10 is found in the external regions 
of C3 and C4. 

Also in the case of clusters our distributions of the ther- 
malised energy flux are steeper than those reported in 
other similar works: we find a t h ~ —4 to —5, while 
a t h ~ —3 to —4 is obtained by Pfrommer et al.(2007), 
where the Lagrangian SPH code Gadget-2 was em- 
ployed. 

The radial profile of the ratio fcn/fth i n our clus- 
ters is reported in Fig[25] Here we show the results in 
the case of both the KJ02 (left panel) and KJ07 (right 
panel) models. Inside the virial radius we do not find 
any relevant difference between our clusters. This is 
because, independently of the dynamical status of the 
clusters, the bulk of the energy is dissipated in the 
CR acceleration and in the thermal component at rel- 
atively weak shocks. The maximum value of fcn/fth 
is found at distance ^ R v i r from the cluster center, 
fcn/fth ~ 0.5 and 0.3 using the KJ02 and KJ07 
model, respectively. 



8.1 Results 

We simulated a large cosmic volume, (145Mpc) 3 m 
(103Mpc/h) 3 , with a fixed grid resolution of 125 kpc. 
Additional simulations were designed and used to in- 
vestigate the effect of spatial resolutions and of the erg 
parameter (see Appendix). 

In the following we summarize the main results ob- 
tained : 

• Re-ionization: in Sect l6.3l we have shown that a 
correct treatment of the re-ionization is crucial to have 
a viable description of shocks with a post processing 
procedure. In particular using a constant temperature 
floor in the simulated data at a fixed redshift may cause 
a somewhat artificial flattening (or pile up) of the dif- 
ferential distribution of the shocks with Mach number. 
This is because the velocity of the sound is artificially 
boosted up in low density and cold regions. To over- 
come this point we derive formulas giving the typical 
temperature of the gas as a function of the local den- 
sity by fitting data obtained from simulations which in- 
clude a specific modeling of the re-ionization in run- 
time. These formulas are found to be consistent with 
Katz et al.(1996) and Valageas et al.(2002) and can be 
used to model the temperature background of adiabatic 
simulations in a post processing procedure. In Section 
4 (Fig© we have shown that our post processing pro- 
cedure is indeed consistent with the results from simu- 
lations with run-time re-ionization. 

• Met hod s to derive the Mach number of shocks: in 
Sect |5.2| and |5.3| we have discussed in some detail two 
different methods to catch shocks in simulated data and 
to estimate their Mach number. These methods are the 
temperature jump (TJ) and the velocity jump (VJ) and 

6 See also Skillman et al.2008, which appeared after this paper was sub- 
mitted. 
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Figure 24. Distribution of thermalised fluxes for the galaxy clusters presented in the text. Distribution are normalized for the volume of the most massive one, 
and are taken from spheres of radius = 2R v i r (left)and R v i r (right) around each cluster. 
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Figure 25. Volume averaged profiles of fen/ fth f° r the four galaxy clusters. Left panel is for the KJ02 model and right panel is for the KJ07 model. 



rely with the jumps in temperature and velocity across 
shocked cells, respectively. 

The shock discontinuity is typically spread over a 
few cells and the risk in measuring the Mach number 
of shocks trough cell-to-cell velocity or temperature 
jumps (or jumps across a few cells) is to underestimate 
the Mach number of shocks. To study this point in the 
context of our numerical simulations we perform sev- 
eral shock-tube tests with ENZO and calculate maps of 
shocks in our simulations under different approaches 
(Sect. 6. 1). We conclude that shocks in our simulations 
are best characterized from velocity (VJ) or tempera- 
ture (TJ) jumps measured across three cells. 

Both the VJ and TJ schemes use ideal conditions 
across non shocked cells (i.e. no velocity and no tem- 
perature gradients) and this may cause major uncer- 
tainties in the characterization of the shocks in a post 



processing procedure (Sect l6.2.~fl - I6.2.21 ). This is be- 
cause the velocity field and temperature distributions 
in the cosmological data sets are very complex and 
the passage of a shock establishes thermodynamical 
gradients which are superimposed to already existing 
ones. We discuss the strength of the uncertainties on 
the value of the Mach number from the two schemes 
by means of Monte Carlo extractions of temperature 
and velocity variations across non shocked cells in our 
data sets, and show that the VJ method may be more 
reliable, at least in the case of weak shocks and espe- 
cially in low density environments (Figs. 8-9). 
Besides these uncertainties we find that the two meth- 
ods yield statistically similar Mach number distribu- 
tions of shocked cells in our simulations (Fig. 12) sug- 
gesting that the characterization of shocks in our sim- 
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ulations is fairly solid; in Sect. 7 we adopt here the VJ 
method. 

• Morphology of LS shocks: in Sect l7.ll we dis- 
cuss the morphology of the shock-patterns detected 
in our simulated data sets. About 15 per cent of the 
cells in our simulations are found to host shocks at the 
present epoch, and this fraction slightly decreases with 
look back time for post-reionization epochs. In qualita- 
tive agreement with previous studies (Ryu et al. 2003; 
Pfrommer et al. 2006) we find that these shocked cells 
form spectacular and complex patterns associated with 
the cosmic web. Filamentary or sheet-like shocks are 
found outside the virial regions of clusters and around 
filaments, while more regular spherical structures sur- 
round galaxy clusters. 

• Number Distributions of LS shocks: we study 
the number distribution of shocked cells as a func- 
tion of the Mach number of the shocks. An impor- 
tant point here is that thanks to the Eulerian scheme 
of the ENZO code we were able to follow the hydro- 
dynamics of the LS shocks also in very low density 
regions, whose exploration is challenging for present 
Lagrangian schemes. 

We find that the bulk of cosmological shocks is es- 
sentially made by weak M ^ 2 shocks and that the 
number distribution of shocks can be grossly described 
by a steep power law MdN/dM oc M a . When con- 
sidering the Mach number distribution of shocked cells 
in the total simulated volume we find an overall steep 
distribution a « —1.6 which is dominated by the 
contribution from voids and filaments. This distribu- 
tion steepens with increasing the cosmic over-density 
and becomes very steep (a « —3 to —4) in the case 
of galaxy clusters demonstrating the increasing rarity 
of strong shocks in these denser (and hotter) regions, 
where basically M < 3. 

• Energy dissipated at LS shocks : the energy dis- 
sipation at LS shocks is the main focus of the pre- 
vious studies on this topic (e.g., Miniati et al.2001; 
Keshet et al.2003; Ryu et al.2003; Pfrommer et al. 
2006; Kang et al.2007; Pfrommer et al.2007). Follow- 
ing Ryu et al.(2003) we calculate the energy rate dis- 
sipated in the form of thermal energy at shock s by 
means of hydrodynamical jump conditions (Eqjl lj). 
In agreement with these previous studies we find that 
about ~ 4 • 10 47 erg/ s are dissipated at shocks in a 
(103Mpc/h) 3 volume in the simulations at the present 
epoch. The bulk of the energy in our simulations is dis- 
sipated in galaxy clusters which indeed contribute to 
about 75 per cent of the total energy dissipation (about 
80 per cent if we also include the contribute from clus- 
ter outskirts), while filaments contribute to a 15 per 
cent of the total energy dissipation. We calculate the 
distribution of the energy flux dissipated at LS shocks 
with shock-Mach number. When considering shocked 
cells in the total volume we find that the distribution is 
steep, a th -2.7 (f th (M)M oc M ath ) and peaks at 
M « 2. Although in qualitative agreement with pre- 
vious studies, our distribution is steeper than that ob- 
tained by Ryu et al.(2003) that also used cosmological 
simulations based on a Eulerian scheme. This differ- 



ence is mostly due to our more complex treatment of 
the re-ionization background and to the use of the VJ 
scheme to measure the Mach number of shocks. 

Following Ryu et al.(2003) we calculate the effi- 
ciency of the injection of CR at LS shocks according 
to Kang & Jones (2002). We obtained Mach number 
distributions of the energy flux dissipated into CR ac- 
celeration in our simulated data sets. Our results are in 
line with previous findings, although our distributions 
are steeper than those in Ryu et al. and slightly steeper 
than those in Pfrommer et al.(2006 & 2007). 

In agreement with Pfrommer et al. we find that the 
bulk of the energy dissipated in the form of CR at 
shocks is shared between clusters and filaments so that 
CR-acceleration happens in regions broader than those 
where thermal energy is dissipated at shocks and the 
ratio between CR energy and thermal energy increases 
in lower density regions (FigfTT]). When considering all 
the shocked cells in our simulations we find that the ra- 
tio between the energy dissipated in the form of CR- 
acceleration and of thermal energy at present epoch 
is fca/fth ~ 0.2 and this ratio becomes smaller in 
galaxy clusters. These fractions are a factor «2 smaller 
than those found by Ryu et al. and this is likely due to 
the steeper distributions that we obtained in our analy- 
sis. 

• Galaxy Clusters: in Sect j7.4| we briefly discuss 
the case of shocks propagating in galaxy clusters. We 
find very steep distributions of the number of shocks 
and of the energy dissipated at shocks as a function 
of the shock-Mach number. The typical Mach number 
within the virial radius is M 1.5, in nice agreement 
with semi-analytical studies which are indeed appro- 
priate for virialised systems (Gabici & Blasi 2003). At 
larger distance from the cluster center stronger shocks 
are found and their presence is strongly correlated with 
the dynamical status of the clusters. Remarkably the 
rarity of moderate-strong shocks in the cluster central 
regions (within rj Mpc distance from cluster center) 
makes the ratio fcii/fth very small, expecially when 
the Kang & Jones (2007) model for the injection of CR 
at shocks is adopted (Fig. 25). 

8.2 Comment on the injection of CR 

The use of numerical simulations is mandatory to un- 
derstand the properties of LS shocks and to investigate 
their role in the injection of CR in LSS. Although we 
use a different approach with respect to other studies, 
our findings for the energy dissipated in the form of 
CR at these shocks are grossly consistent with previ- 
ous studies (Ryu et al. 2003; Pfrommer et al. 2006). 

However, the astrophysical problem is extremely 
complex and several hidden ingredients in the adopted 
procedures are potentiall y sou rces of large uncertain- 
ties. As discussed in Sect j7.3l the efficiency of CR ac- 
celeration at shocks is investigated following several 
approaches. In the present paper we have used the ac- 
celeration efficiency resulting from numerical calcula- 
tions of modified shocks (following Ryu et al. 2003 
and references therein). On the other hand Pfrommer et 
al.(2006) use a linear theory with the efficiency mod- 
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ified to account for saturation effects at large values 
of the Mach numbers (actually to limit the CR effi- 
ciency at m 50 per cent). These two approaches are 
formally radically different, but nevertheless they pro- 
vide an overall estimate of the CR injection efficiency 
which is not dramatically different in the case of typ- 
ical shocks with M » 2 - 4. The main hidden ingre- 
dient in the efficiency of CR acceleration comes from 
the commonly adopted thermal leakage injection sce- 
nario which essentially adopts as minimum momentum 
of the particles that take part in the acceleration pro- 
cess, Pinj, a multiple of the momentum of the thermal 
particles, pi n j = XiPth- The choice of Xi is a guess, 
since this depends on physical details which are still 
poorly known (e.g., Blasi 2004). In Ryu et al.(2003) 
(and thus in our paper) the fraction of protons injected 
into the CR population at shocks is taken of the or- 
der of « 10~ 3 which is not far from (even if larger 
than) the resulting efficiency from the assumption of 
p in j 3.5pth adopted in Pfrommer et al.(2006). Al- 
though this parameter is somewhat constrained by the 
theory (e.g., Malkov 1998), it should be stressed that 
having a slightly different value of xi (e.g. £j=3.8 in- 
stead of 3.5) would have the net effect to reduce the 
acceleration efficiency by nearly one order of magni- 
tude. 



8.3 Constraints from observations 

Theoretical arguments suggest that the bulk of CR in 
galaxy clusters should be in the form of supra-thermal 
protons (e.g. Blasi, Gabici, Brunetti 2007 for a re- 
cent review). Gamma ray observations of a few nearby 
galaxy clusters limit the energy density of CR protons 
in these clusters to « 20 per cent of the thermal energy 
(Pfrommer & Ensslin 2004; Reimer 2004). Although a 
detailed comparison with these limits would require to 
follow the advection and accumulation process of CR 
in clusters during the simulations (e.g. Miniati 2003; 
Pfrommer et al. 2006), the level of injection rate o f CR 
that we found in clusters and cluster outskirts (Sect l7.4l ) 
allows us to reasonably conclude that our results are 
barely consistent with these limits. 

On the other hand more stringent upper limits can 
be obtained from present radio observations. The bulk 
of galaxy clusters does not show evidence of extended 
Mpc-scale synchrotron radio emission and this can be 
used to constrain the population of secondary electrons 
and thus that of the primary CR protons from which 
these secondaries would be injected (Brunetti et al. 
2007). These limits are very stringent and actually rep- 
resent a challenge for simulations: in the case that the 
ICM is magnetized at « /j,G level (consistent with Ro- 
tation Measures, e.g. Govoni & Feretti 2004) the en- 
ergy of CR should be at ^ few percent of the thermal 
energy (when the spectrum of CR is fixed at that ex- 
pected from simulations, i.e. s « 2 — 2.5, N(p) oc p~ s 
for M > 3). 

A comparison between our simulations and present 
limits clearly requires a more detailed study. However, 
a simple estimate of the spectral shape of CR injected 



in our simulations (Figl26l) suggests that the bulk of the 
CR energy in clusters and cluster outskirts is associ- 
ated with CR populations with relatively steep spectra, 
in which case the limits from radio observations be- 
come less stringent (Brunetti et al. 2007). Obviously, 
an alternative possibility to reconcile radio data and 
simulations would be that relatively flat components 
of CR store a non negligible fraction of the energy of 
the ICM, and that the magnetic field strengths in the 
ICM are much smaller than that claimed from the anal- 
ysis of present RM data; future gamma ray observa- 
tions with the FERMI Gamma Ray Telescope and with 
Cherenkov arrays will clarify this possibility. 
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numerical resolutions (i.e. AD125, AD250, AD500 and ADS00). 



Volk H. J., Aharonian F. A., Breitschwerdt D., 1996, 
SSRv, 75, 279 

Woodward P., Colella P., 1984, JCoPh, 54, 115 



APPENDIX 



APPENDIX A: THE EFFECT OF SPATIAL RESOLUTION. 

We investigated the effect of resolution on the proper- 
ties of detected shocks by re-simulating the same ini- 
tial conditions and cosmic volume of the AD 125 sim- 
ulations at resolutions of 800/cpc, 500kpc and 250kpc. 

Even if most of the graphs and statistics presented 
in the paper are done by using n = 1 for the shock 
detecting scheme (see Sec.5.3) and thus assuming that 
the best reconstruction of the shock discontinuity is 
achieved by considering a jump of 2 cells between 
pre-shock and post-shock, here we prefer to keep this 
jump smaller (i.e. n = 0). This is in order minimize 
any confusion coming from the fact that in poorly re- 
solved runs shocks have sizes of typical cluster halos 
(i.e. for n = 1 in the AD800 one would reconstruct 
shocks across 1.6Mpc). This is fair enough to recon- 
struct the trend with resolution within our simulations, 
and the comparison to the n = 1 case can be recovered 
in Fig. 14. 

First we find that at all these resolutions Eq. [Qpro- 
vides a good fit to the density-temperature distribu- 
tions obtained with run time re-ionization and thus we 
use this relation to model the re-ionization in our post 
processing approach at all these resolutions. We then 
analyze the outputs at z = and derive statistical prop- 
erties of shocks in the simulated volumes, following 
the procedures given in the previous Sections. 
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Figure A2. Distribution of the thermalised energy flux in different over- 
density bins, for 4 different numerical resolution. 



Figure Bl. Effect of the variation of the erg parameter on the distribution of 
thermalised flux, in different over-density bins. 



The number distributions of shocked cel ls a s a 
function of their Mach number are given in Fig. lAll for 
the different resolutions. We find that the results con- 
verge at higher resolutions, in particular the shape of 
the distribution and integral number of shocked cells 
obtained with I25kpc and with 250/cpc resolution are 
consistent within 20 percent. A relevant point here 
is that the excess of shocks with high Mach number 
found at low resolution is progressively reduced with 
increasing resolution. 

The case of the ener gy flux dissipated at shocks is 
reported in Figures IA2I This case is somewhat more 
problematic as it depends on the combination of the 
properties of shocks with the local baryon over-density. 
Despite the properties of shocks statistically converge 
with resolution, the over-density in the simulated vol- 
ume increases with resolution and this causes the in- 
crease of the dissipated energy in the data simulated 
at higher resolutions. Anyhow also in this case some 
level of convergence is obtained in line with previous 
studies (Ryu et al. 2003; Pfrommer et al. 2006a). The 
most problematic case is that of galaxy clusters were 
indeed the dissipation of the energy at shocks increases 
by one order of magnitude between lower and higher 
resolution data sets (this still increases by rj 1.5 times 
between the 250 and 125 kpc data sets). 

Despite this slow convergence with resolution, the 
value of the ratio Jcr/ fth is found to not significantly 
change with resolution since the resolution affects the 
two quantities in a similar way. 



APPENDIX B: THE EFFECT OF A VARIATION OF THE 
a 8 PARAMETER. 

The value of the eg parameter (the normalization in 
the power spectrum of primordial over-density fluc- 
tuations) crucially affects the abundance of collapsed 
objects in the universe at a given epoch. This value is 
presently not fully constrained: very recent CMB anal- 
ysis give a relatively small value, erg = 0.74 (Spergel 
et al.2007), with respect to that derived from previous 



CMB data-analysis (Spergel et al.2003) and to the con- 
straints from the observed abundance of galaxy clusters 
(e.g. Evrard et al.2007). In this Appendix we briefly 
discuss the effect of the <jg parameter on the statistical 
properties of the shocks as measured in our simulated 
data-sets (adopting as in the previous Section n = for 
the reconstruction of shocks). We thus re simulated the 
CO250 run with a 8 = 0.74 (S8250) and applied all the 
procedures discussed in the previous Sections to derive 
the properties of the shocks (note that the CO250 and 
S8250 simulations have run-time re-ionization). 

Theoretically the population of shocks in a uni- 
verse with larger <7g is expected to evolve faster as more 
power is associated with the primordial over-density 
fluctuations. Thus, at a fixed redshift, universes with 
larger erg host more evolved structures, which are char- 
acterized by typically higher internal sound speeds at 
higher densities, and low temperatures in low density 
regions. 

The distribution of thermalised energy at shocks 
in the two simulations is given in Fig.Bl. Although 
clearly modifying the value of <jg has some effect on 
the properties of the shocks in the simulations, the net 
result is that, within the presently allowed region of the 
values of the <jg parameter, no strong difference in the 
properties of the shocks is found in simulations with 
different <7g. Indeed, globally we find that the energy 
dissipated at the present time in the S8250 simulation 
is 2 times smaller than that in the CO250 simula- 
tion, and the distribution with Mach number of the dis- 
sipated energy in under-dense regions is found to be- 
come somewhat slightly flatter with decreasing <7g. 



